 ---
title: "Modulation Ongoing Oscillations"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r, include=FALSE}
library(rstudioapi)  #to set working directory automatically to file location
library(matrixStats)
library(ggplot2) #plots
library(ggpubr) #publiation ready plots
library(base)
library(agricolae)
library(sjmisc)
library(psych)
library(irr)
library(rstatix)
library(broom)
library(purrr)
library(lmerTest) #for anova
library(stats) # for anova
library(lme4) #for lmer anova
library(report) #to report lmer automatically
library(emmeans) #pairwise comparisons
library(pwr)
library(simr) #for sample size / power simulation
library(bestNormalize)
library(ggsci) #coloring of plot >10 colors
library(car)
library(pbkrtest) #for kenward-rogers approximation of LMM results
library(sjPlot) #for nice table of LMM (function tab_model())
```


### load data

```{r, include=FALSE}
setwd(dirname(getActiveDocumentContext()$path))       # Set working directory to source file location

theta_modulation <-read.delim("aggregatedHarmonics_newICA_theta.csv", header=TRUE, sep=";")             #import dataset
theta_modulation <-convert_as_factor(theta_modulation, vars=  c("modality", "condition", "electrode"))  #set variables as factors
```

### sort table

```{r}
#selects electrode of interest (defined by Wilcoxon test), if the same for both modalities
electrodeOI <-  as.data.frame(filter(theta_modulation, electrode == "C3"  )) 

#selects electrode of interest (defined by Wilcoxon test), if different for different modalities
electrodeOI2 <-  as.data.frame(filter(theta_modulation, electrode == "CP3" & modality == "thermonociceptive" &  electrode == "CP4" & modality == "vibrotactile" ))

electrodeOI$Subject<-as.factor(rep(sprintf("Sub_%d", 1:25),times=2))  #add factor for #Subjects
electrodeOI <- convert_as_factor(electrodeOI, vars=  c("modality", "condition", "electrode")) #convert to factors
```


#### Removing artifacted subject from dataset
```{r}
electrodeOI<- subset(electrodeOI, Subject!="Sub_17") #removes all the rows containing "Sub_17" (which is Subject 25 in reality)
```

### Outlier Tests

#### Primary outlier check (Least Square )
```{r}
#create LMM for outlier removal
LMM_amplitudeMod_pre <- lmer(amplitude ~ condition + (1|Subject), data=electrodeOI)

# remove outlier
#modulation_clean <-modulation$amplitude[which(modulation$amplitude > Tmin & modulation$amplitude < Tmax)] 

#look at data for outliers
boxplot(electrodeOI$amplitude)

# get mean and Standard deviation
mean = mean(electrodeOI$amplitude)
std = LMM_amplitudeMod_pre@sigma

# get threshold values for outliers
Tmin = mean-(3*std)
Tmax = mean+(3*std)

# find outlier and replace with NA directly in dataset
electrodeOI$amplitude <- replace(electrodeOI$amplitude, c(electrodeOI$amplitude< Tmin | electrodeOI$amplitude > Tmax),NA)
```


#### Assumptions and influential data points
Testing if assumption of normality and linearity is respected and whether there is a data point which overproportionally influences the dataset (Cook's Distance >1)
```{r}
lmm <- lmer(amplitude ~ condition*modality+ (1|Subject), data=electrodeOI) #testing assumptions with Linear Mixed Model 

#creating a QQplot to test the assumption of normality
qqnorm(resid(lmm))
qqline(resid(lmm))

#to check if a datapoint is influental, look at  Cook's distance. If between 0.5<x>1, not necessarly influential. If >1, influential and should be removed
electrodeOI$cook <- cooks.distance(lmm)
plot(electrodeOI$cook, electrodeOI$Subject)


#log transformation to make data more conform with normality, +1 to make amplitude big enough for log to have an effect.
electrodeOI$amp_log <- log10(electrodeOI$amplitude*10+1)
lmm_log <- lmer(amp_log ~ condition*modality+ (1|Subject), data=electrodeOI)
qqnorm(resid(lmm_log))
qqline(resid(lmm_log))

electrodeOI$cook <- cooks.distance(lmm_log)
plot(electrodeOI$cook, electrodeOI$Subject)

shapiro.test(electrodeOI$amp_log)
```

#### Removal of Subject if detected as outlier / not conforming to assumptions
```{r}
#replaces data points of identified Subjects with a "NA"
#electrodeOI$amp_log <- replace(electrodeOI$amp_log,Subject== "Sub_21" | Subject == "Sub_11" ,NA) # example to remove Subject 1 and 5
```


## Linear Mixed Model
```{r}
LMM_amplitudeMod_cook <- lmer(amp_log ~ condition *modality+ (1|Subject), data=electrodeOI,REML = TRUE)

if(requireNamespace("pbkrtest", quietly = TRUE))           #"activating" Kenward Roger approximation
anova(LMM_amplitudeMod_cook, ddf="Kenward-Roger", type =3) #to get LMM test statistic approximated by Kenward Rogers
emmeans(LMM_amplitudeMod_cook, list(pairwise ~ condition| modality), adjust = "tukey" ) #post-hoc pairwise comparison, if necessary
```


### Post hoc power calculation 
```{r}
LMM_amplitudeMod_cook <- lmer(amplitude ~ condition *modality+ (1|Subject), data=electrodeOI)
sim_power<- powerCurve(LMM_amplitudeMod_cook,test=fixed("modality", method="anova"),along="Subject", nsim=100, alpha=0.02, seed=135)
plot(sim_power) #plotting the calculated power for the different sample sizes

sim_power_long <- extend(LMM_amplitudeMod_cook, along = "Subject", n=50) #extending the sample size artificially
sim_power_long<- powerCurve(sim_power_long,test=fixed("modality", method="anova"),along="Subject", nsim=100, alpha=0.02, seed=135) #condition:modality for interaction
plot(sim_power_long) 
```

Dotplot
```{r, echo=FALSE, warning =FALSE}
electrodeOI_plot <-ggerrorplot(electrodeOI, x="modality", y="amplitude",facet.by="condition",combine = TRUE, size= 1, ylim=c(-0.05,0.25),
            add=c( "jitter"), title= FALSE, xlab= FALSE, ylab="theta amplitude [µV]", color="modality", palette = c("red", "blue"),
            panel.labs.background=list(color="white", fill="white"), panel.labs.font=list(face="bold", size=12))+
    theme(panel.border = element_blank(), axis.line = element_line()) +
    geom_bracket( xmin = c("thermonociceptive"), xmax = c("vibrotactile"), y.position = c(0.2), label = c("*"), tip.length = 0.02 ) #dotplot of data, setting statistical significance manually, including confidence intervals

ggadd(electrodeOI_plot,add ="mean_ci", color="black", shape= 17, size= 0.8) #add group mean and confidence interval
```
